#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <string>
#include <cstring>
#include <vector>
#include <stdint.h>

#include "galois-2.h"
#include "jerasure-2.h"
#include "reed_sol-2.h"

static int prim08 = -1;
static int mask08_1 = -1;
static int mask08_2 = -1;

static int prim16 = -1;
static int mask16_1 = -1;
static int mask16_2 = -1;

static int prim32 = -1;
#define rgw32_mask(v) ((v) & 0x80000000)

pthread_mutex_t prim_lock;

class JER_Matrix *reed_sol_vandermonde_coding_matrix(int k, int m, int w)
{
  int tmp;
  int i, j, index;
  JER_Matrix *jm, *vand;

	/* strip the last m rows off of the big vandermonde matrix */
  vand = reed_sol_big_vandermonde_distribution_matrix(k+m, k, w);
  if (vand == NULL) return NULL;
	try {
		jm = new JER_Matrix(m, k, w);
	} catch (bad_alloc &e) {
		cerr << e.what() << endl;
		delete vand;
		return NULL;
	}

  for (i = 0; i < m; i++) {
    for (j = 0; j < k; j++) {
      jm->Elts[i][j] = vand->Elts[i+k][j];
    }
  }
  delete vand;
  return jm;
}
class JER_Matrix *reed_sol_extended_vandermonde_matrix(int rows, int cols, int w)
{
  JER_Matrix *jm;
  int i, j, k;

  if(w < 30 && (1 << w) < rows) return NULL;
  if(w < 30 && (1 << w) < cols) return NULL;

	try {
		jm = new JER_Matrix(rows, cols, w);
	} catch (bad_alloc &e) {
		cerr << e.what() << endl;
    return NULL;
  }

	/* Initialize the first and last rows */
  jm->Elts[0][0] = 1;
  for(j = 1; j < cols; j++) jm->Elts[0][j] = 0;
  if(rows == 1) return jm;

  for(j = 0; j < cols-1; j++) jm->Elts[rows-1][j] = 0;
  jm->Elts[rows-1][j] = 1;
  if(rows == 2) return jm;

	/* Fill in eac row with the multiples of the row index */
  for(i = 1; i < rows-1; i++) {
    k = 1;
    for(j=0; j<cols; j++) {
      jm->Elts[i][j] = k;
      k = galois_single_multiply(k, i, w);
    }
  }
  return jm;
}
class JER_Matrix *reed_sol_big_vandermonde_distribution_matrix(int rows, int cols, int w)
{
  JER_Matrix *jm;
  int i, j, k;
  int sindex, srindex, siindex, tmp;
  vector <unsigned int> vtmp;

  if(cols >= rows) return NULL;

  jm = reed_sol_extended_vandermonde_matrix(rows, cols, w);
  if(jm == NULL) return NULL;

  for(i = 1; i < cols; i++) {

    for(j = i; j < rows && jm->Elts[j][i] == 0; j++);
    if (j >= rows) {
      fprintf(stderr, "reed_sol_big_vandermonde_distribution_matrix(%d,%d,%d) - couldn't make matrix\n",
          rows, cols, w);
      exit(1);
    }

    /* If necessary, swap rows */
    if(j != i) {
      vtmp = jm->Elts[i];
      jm->Elts[i] = jm->Elts[j];
      jm->Elts[j] = vtmp;
    }

    /* If element i,i is not equal to 1, multiply the column by 1/i */
    if(jm->Elts[i][i] != 1) {
      tmp = galois_single_divide(1, jm->Elts[i][i], w);
      for(j = 0; j < rows; j++) 
        jm->Elts[j][i] = galois_single_multiply(tmp, jm->Elts[j][i], w);
    }

    /* Now, for each element in row i that is not in column 1, you need
       to make it zero.  Suppose that this is column j, and the element
       at i,j = e.  Then you want to replace all of column j with 
       (col-j + col-i*e).   Note, that in row i, col-i = 1 and col-j = e.
       So (e + 1e) = 0, which is indeed what we want. */

    for(j = 0; j < cols; j++) {
      tmp = jm->Elts[i][j];
      if(j != i && tmp != 0) {
        for(k = 0; k < rows; k++) {
          jm->Elts[k][j] = jm->Elts[k][j] ^ galois_single_multiply(tmp, jm->Elts[k][i], w);
        }
      }
    }
  }

  /* We desire to have row k be all ones.  To do that, multiply
     the entire column j by 1/Elts[k][j].  Then row j by 1/Elts[j][j]. */

  for (j = 0; j < cols; j++) {
    tmp = jm->Elts[cols][j];
    if (tmp != 1) {
      tmp = galois_single_divide(1, tmp, w);
      for (i = cols; i < rows; i++) {
        jm->Elts[i][j] = galois_single_multiply(tmp, jm->Elts[i][j], w);
      }
    }
  }

  /* Finally, we'd like the first column of each row to be all ones.  To
     do that, we multiply the row by the inverse of the first element. */

  for (i = cols+1; i < rows; i++) {
    tmp = jm->Elts[i][0];
    if (tmp != 1) {
      tmp = galois_single_divide(1, tmp, w);
      for (j = 0; j < cols; j++) jm->Elts[i][j] = galois_single_multiply(jm->Elts[i][j], tmp, w);
    }
  }
  return jm;
}

class JER_Matrix *reed_sol_r6_coding_matrix(int k, int w) {
  JER_Matrix *jm;
  int i, tmp;

  if (w != 8 && w != 16 && w != 32) return NULL;
  
	try {
		jm = new JER_Matrix(2, k, w);
	} catch (bad_alloc &e) {
		cerr << e.what() << endl;
		return NULL;
	}

	/* first row is all ones */
  for (i = 0; i < k; i++) jm->Elts[0][i] = 1;
	/* indices of the row are mutliples of two */
  tmp = 1;
  for (i = 0; i < k; i++) {
    jm->Elts[1][i] = tmp;
    tmp = galois_single_multiply(tmp, 2, w);
  }

  return jm;
}

int reed_sol_r6_encode(class JER_Slices *slices)
{
  int i, j, k;
  int pstarted;
  int sindex;
  unsigned char *dptr, *pptr, *bdptr, *bpptr;

  /* First, put the XOR into coding region 0 of all slices */
  bdptr = slices->ptrs[0];
  bpptr = slices->ptrs[slices->K];
  memcpy(bpptr, bdptr, slices->size);

  for (i = 1; i < slices->K; i++) galois_region_xor(bpptr, slices->ptrs[i], bpptr, slices->size);

  /* Next, put the sum of (2^j)*Dj into coding region 1 */
  bpptr = slices->ptrs[slices->K+1];
  for (sindex = 0; sindex < slices->size; sindex+=slices->PS) {
    pptr = bpptr + sindex;
    dptr = slices->ptrs[slices->K-1] + sindex;
    memcpy(pptr, dptr, slices->PS);

    for (i = slices->K-2; i >= 0; i--) {
      dptr = slices->ptrs[i] + sindex;
      switch(slices->W) {
        case 8 : reed_sol_galois_w08_region_multby_2(pptr, slices->PS); break;
        case 16: reed_sol_galois_w16_region_multby_2(pptr, slices->PS); break;
        case 32: reed_sol_galois_w32_region_multby_2(pptr, slices->PS); break;
        default: return 0;
      }

      galois_region_xor(pptr, dptr, pptr, slices->PS);
    }
  }
  return 1;
}

void reed_sol_galois_w08_region_multby_2(unsigned char *region, int nbytes)
{
  uint32_t *l1;
  uint32_t *ltop;
  unsigned char *ctop;
  uint32_t tmp, tmp2;

	pthread_mutex_lock(&prim_lock);
  if (prim08 == -1) {
    tmp = galois_single_multiply((1<<7), 2, 8);
    prim08 = 0;
    while (tmp != 0) {
      prim08 |= tmp;
      tmp = (tmp << 8);
    }
    tmp = (1 << 8) - 2;
    mask08_1 = 0;
    while (tmp != 0) {
      mask08_1 |= tmp;
      tmp = (tmp << 8);
    }
    tmp = (1 << 7);
    mask08_2 = 0;
    while (tmp != 0) {
      mask08_2 |= tmp;
      tmp = (tmp << 8);
    }
  }
	pthread_mutex_unlock(&prim_lock);

  ctop = region + nbytes;
  ltop = (uint32_t *) ctop;
  l1 = (uint32_t *) region;

  while (l1 < ltop) {
    tmp = ((*l1) << 1) & mask08_1;
    tmp2 = (*l1) & mask08_2;
    tmp2 = ((tmp2 << 1) - (tmp2 >> 7));
    *l1 = (tmp ^ (tmp2 & prim08));
    l1++;
  }
}
void reed_sol_galois_w16_region_multby_2(unsigned char *region, int nbytes)
{
  uint32_t *l1;
  uint32_t *ltop;
  unsigned char *ctop;
  uint32_t tmp, tmp2;

  pthread_mutex_lock(&prim_lock);
  if (prim16 == -1) {
    tmp = galois_single_multiply((1 << 15), 2, 16);
    prim16 = 0;
    while (tmp != 0) {
      prim16 |= tmp;
      tmp = (tmp << 16);
    }
    tmp = (1 << 16) - 2;
    mask16_1 = 0;
    while (tmp != 0) {
      mask16_1 |= tmp;
      tmp = (tmp << 16);
    }
    tmp = (1 << 15);
    mask16_2 = 0;
    while (tmp != 0) {
      mask16_2 |= tmp;
      tmp = (tmp << 16);
    }
  }
	pthread_mutex_unlock(&prim_lock);

  ctop = region + nbytes;
  ltop = (uint32_t *) ctop;
  l1 = (uint32_t *) region;

  while (l1 < ltop) {
    tmp = ((*l1) << 1) & mask16_1;
    tmp2 = (*l1) & mask16_2;
    tmp2 = ((tmp2 << 1) - (tmp2 >> 15));
    *l1 = (tmp ^ (tmp2 & prim16));
    l1++;
  }

}
void reed_sol_galois_w32_region_multby_2(unsigned char *region, int nbytes)
{
  uint32_t *l1;
  uint32_t *ltop;
  unsigned char *ctop;

	pthread_mutex_lock(&prim_lock);
  if (prim32 == -1) prim32 = galois_single_multiply((1 << 31), 2, 32);
	pthread_mutex_unlock(&prim_lock);

  ctop = region + nbytes;
  ltop = (uint32_t *) ctop;
  l1 = (uint32_t *) region;

  while (l1 < ltop) {
    *l1 = ((*l1) << 1) ^ ((*l1 & 0x80000000) ? prim32 : 0);
    l1++;
  }
}
